Evolution of the macroscopically entangled states in optical lattices 
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We consider dynamics of boson condensates in finite optical lattices under a slow external pertur- 
bation which brings the system to the unstable equilibrium. It is shown that quantum fluctuations 
drive the condensate into the maximally entangled state. We argue that the truncated Wigner 
approximation being a natural generalization of the Gross-Pitaevskii classical equations of motion 
is adequate to correctly describe the time evolution including both collapse and revival of the con- 
densate. 
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I. INTRODUCTION 



Recent advances in experimental realization of Bose- 
Einstein condensates (BECs) in optical latticesA* 2 ^ make 
this field particularly interesting for theoretical analy- 
sis. One of the most striking features about these con- 
densates is the possibility to observe directly effects of 
quantum fluctuations at zero temperature. For exam- 
ple, as was predicted theoretically^ and shown experi- 
mentally, the zero point motion can drive the system 
from the superfiuid to the Mott insulating state. The 
other direct manifestation of the quantum effects was re- 
ported in Ref.Q, where it has been shown that bosons 
can live in the superposition of number states even at 
the absence of tunnelling. On the other hand, in the su- 
perfiuid regime the quantum fluctuations are suppressed 
and either classical Gross-Pitaevskii (GP) or Bogoliubov 
approach is often adequate for the description of both 
static and dynamic properties of the condensates (see 
e.g. Refs. jalflj). However, there is an interesting pos- 



sibility, wherein the system is superfiuid but neither of 
these approaches is good. Suppose that the initially sta- 
ble condensate is driven to the regime of instability. This 
can be achieved either by applying a certain phaseshift 
to the condensate with repulsive interactions^^ or by 
switching the sign of the interaction to the negative value 
using Fcshbach resonanceiSiii . The main difficulty arises 
because near the instability all the fluctuations including 
quantum exponentially diverge and cannot be treated as 
a small perturbation. To be more specific, suppose that 
for time t < the periodic system of condensates in a 
lattice was in a superfiuid ground state, i.e. the interac- 
tion was relatively weak. Then a phase imprint, i.e. a 
certain phase difference between the adjacent wells, was 
imposed. Experimentally it can be achieved by e.g. ap- 
plying a short (compared to a single tunneling time) pulse 
of external field to the system. A case of special inter- 
est will be when there is a relative it phase shift between 
neighboring wellsa*. For the two wells with equal number 
of bosons and relatively small interaction , this state is 
metastable^i 2 (this is also the case for even number of 
wells and periodic boundary conditions). However, if the 
interaction increases becomes larger than a critical value, 
this equilibrium becomes unstable and the bosons spon- 
taneously form a "dipole" state2ii2iiii4 in which most 



of them occupy one of the two wells. Upon accounting 
for quantum fluctuations in a system with a finite num- 
ber of bosons, the state obtained is a superposition of 
the two dipole states so that the inversion symmetry is 
preserved. Clearly in the case of infinite number of wells 
translational symmetry is always broken. For example in 
Ref. [l^ a similar instability but for the case of a Mott 
insulator in a strong electric field was shown to drive the 
system into a dipole state. 

Related to this instability is a very interesting possibil- 
ity of forming a Schrodinger cat states*. If the interaction 
slowly increases in the 7r state, then as we just mentioned, 
at certain point the system becomes unstable. Classically 
the bosons will remain in this unstable state forever un- 
less there is some noise present either dynamical or in the 
initial conditions. As we will show below such a noise will 
drive all the bosons into one spontaneously chosen well. 
However, apart from classical fluctuations, which are al- 
ways there but relatively weak in the condensates, there 
is also a quantum zero point motion, which comes from 
the uncertainty relation between the number of bosons 
and their phase, so that the state where both are defined 
is simply impossible. This quantum noise will also cause 
the classical trajectories to move apart from the unsta- 
ble equilibrium. However, as we mentioned above, the 
quantum fluctuations do not break translational invari- 
ance so the resulting state must be macroscopically en- 
tangled. Let us give a simple analogy with a ball laying 
on the top of the hill. Without fluctuations it will re- 
main there forever. However, because of the uncertainty 
principle this ball will move down along different classi- 
cal paths. The quantum effects will be manifested only 
in certain phase relations between these paths but will 
not affect the motion itself. This analogy suggests that a 
good way to describe these situations is to take into ac- 
count fluctuations yielding some probability distribution 
of the number and phase at t = and evolve the fields 
according to the classical equations of motion. In the lit- 
erature this approach is known as the truncated Wigner 
approximation (TWA) 17 i 18 i 19 i 20 i 21 i 22 . In the Appendix |X| 
we will show that this approach naturally arises in the 
path integral derivation of the evolution equations, and 
in Sec. IHII we will numerically test the results on the ex- 
actly solvable model of two condensates. We show how 
to go beyond TWA in a separate publication 27 . 

Let us also briefly mention a few other alternatives to 
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generalize the classical dynamics of the BECs existing in 
literature. One of them is a generalization of GP equa- 
tions by adding the interaction of the superfluid compo- 
nent of the condensate with excited bosons2i2i. This 
method was successfully applied for the description of 
the condensate evaporation after a sudden change in the 
scattering length. It is hard to see though what is the 
range of applicability of the resulting equations and how 
crucial are the assumptions made about the dynamics of 
incoherent excitations. Recently, there was developed an- 
other class of methods based on exact stochastic dynam- 
icsi&25iSi. These ideas look very promising, but so far 
the theory is applicable to a particular class of two-body 
interactions and at least to author it is not completely 
clear how to generalize it. 

Throughout this paper we will explicitly consider a one 
dimensional array of coupled condensates with. However, 
the results are quite general and should not depend on 
dimensionality. 

The standard Bose-Hubbard hamiltonian we are going 
to employ reads: 

n = Y^ -J{a]a j+1 + a) +1 aj) + ^ <'>,(<'>; - 1). (1) 
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Here a,j is the canonical Bose annihilation operator on 
sites of the optical lattice (wells) labelled by an integer j, 
J is the tunneling amplitude between neighboring lattice 
sites, U > is the repulsive interaction energy between 
bosons in the same well. Another important parameter 
in the problem is the mean number of bosons per lattice 
site - N. Throughout this paper we consider the case of 
large N, since it corresponds to the nearly classical limit. 
A dimensionless measure of the strength of interactions 
between the bosons is the coupling^ 



Hereafter, except otherwise specified, we set ft = 1 and 
J = 1 so all the energies are given in units of J and time 
has units of ft,/ J. As we noted beforei^, A ~ 1 corre- 
sponds to the crossover from weakly to strongly interact- 
ing superfluid and A ~ ./V 2 corresponds to the quantum 
phase transition to the Mott insulating phase. We will 
be interested only in the superfluid regime and assume 
that A <C N 2 . Note that the hamiltonian (|TJ clearly 
has a time reversal invariance. Besides the equations of 
motion: 

i^- = -[H, ai ] (3) 

are invariant under the transformation: t — > — t, A — > —A, 
and a,j — ► (— l) J dj- So in the absence of energy relax- 
ation, which would break the time reversal symmetry, a 
7r phase-shift between neighboring sites is equivalent to 
the change of the sign of the interaction from repulsive to 
attractive. This equivalence is very useful for qualitative 
understanding of the resulting instabilities. 



II. SEMICLASSICAL EVOLUTION OF THE 
PHASE MODULATED STATE 

As we discussed in detail in Ref. 0]> in the classical 
description of the two coupled condensates the effective 
motion of the number difference is equivalent to that of 
a classical particle with a unit mass in the effective po- 
tential: 




. (4) 

where the "coordinate" n represents the number differ- 
ence between the left and the right sites, 9 is the relative 
phase; no = n(t = 0) and 9o = 0{t — 0) are the initial 
conditions. In particular, if n — and 6> = ir then 

C/ e// = 2n 2 (l-A) + ^. (5) 

Clearly the equilibrium n — becomes unstable if A > 
A c = 1. 

In a more general case of multiple wells a similar anal- 
ysis can be done. Because there are now many degrees 
of freedom a simple representation of the motion using 
the effective potential becomes impossible. Instead let us 
return to the GP version of the equations |JSJ : 

where ipj(t) is the semiclassical field corresponding to 
the expectation value of the operator dj(t). Clearly the 
7r modulated state is a stationary solution for any inter- 
action: 

Ut) = (-l)V e «. (7) 
The unimportant global phase <3(t) is given by: 

9(f) = -2t- f A(r)dr. (8) 
Jo 

Similarly to the two well case this state becomes unstable 
when the interaction exceeds a certain critical valu ofi 12 : 

A c = 2sin 2 ^:, (9) 

where M is the number of the lattice sites in the ar- 
ray. The origin of the instability becomes intuitively 
clear if we use a dynamical symmetry mentioned above: 
A — > —A, tpj — > (—ly-ip*. So the strong repulsive interac- 
tion for the 7r state is equivalent to the strong attractive 
interaction for the symmetric state. The instability for 
the attractive interaction is naturally expectecUfi*ii. To 
get more quantitative results we consider a time evolu- 
tion of fluctuations around the it state: 

^•(t) = (-l)'V e W(l + &(t) +%(*)), (10) 
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with £j and r\j being the small real deviations from the 
exact 7r solution found above. Substituting (fTUf> into © 
and linearizing the resulting equations we obtain: 



-dt = ^ + 



2 Vj 



2A(t)& 



(11) 
(12) 



In the Fourier space this system is equivalent to a set of 
decoupled second order differential equations: 



dt 2 



-16sin 4 |e g + 8A(t)sin 2 |^, (13) 



to study the nonlinear regime of GP equations. More 
specifically the relation 



l£o 9 |e s 



1 



(20) 



defines the boundary between the regimes of small and 
large fluctuations. Using the estimate of £o<j (see the next 
section for the details) 



1 



N 



(21) 



we derive that the instability for a given momentum 
mode will evolve into the nonlinear regime given that 



and 



dH q 



4 sin 2 f dt 



(14) 



In the case of adiabatically changing interaction we can 
write: 



(15) 



and neglect by the second derivative of <f> q . Here £oq is 
the initial amplitude of fluctuations. Substituting l|15|) 
into (|12|) we get: 



dt 



9 =±4sin 2 f< 



2sin 2 f 



(16) 



r 27rsin 2 | 

8 < 2 -. 

~ In TV 



(22) 



This tells us that the interaction should indeed change 
slowly in time (at least near the onset of instability) and 
so justifies the adiabatic limit we used. The lowest energy 
mode corresponds to the momentum q — 2tt/M so the 
lower boundary for 8 becomes: 



Si = 



2lT . 2 7T 

irWV Sm M' 



and the upper boundary is 



2tt 



(23) 



(24) 



For simplicity we assume that the interaction increases 
in time as: 



A(i) = 



A 



l -St' 



(17) 



where 8 is the parameter of adiabaticity and Ao is the 
initial interaction, which we assume to be small . This 
type of A(i) dependence, in fact, corresponds to the tun- 
nelling exponentially decreasing in time (J(t) — Joe~ A *) 
with t — > t/J(t). It is straightforward to solve (|16|l ex- 
plicitly and the result is: 



If 8 < 8i then all the modes have enough time to get into 
the nonlinear regime. On the contrary if 6 > 82 then the 
fluctuations around GP state will remain small. In the 
intermediate regime 8\ < 8 < 82 some of the momen- 
tum modes will exhibit small fluctuations and some will 
become strongly enhanced. 



III. QUANTUM FLUCTUATIONS 
A. Truncated Wigner approximation 



<MA) = ±- 



2 sin 



2 q 



2A 

~ — 



2A 



In- 



1 - 



2sin a a 



1-Jl- 



A 



An 



In- 



1 - 



l-Jl 



Ao 

2 sin^ I -I 



(18) 



Assuming that Ao < 2 sin 2 q/2 we see that in the limit 
A — > 00 the imaginary part of phase 4> q goes to: 



30 9 (oo) 



, 2vr . 2 q 
±— sin -. 

8 2 



(19) 



So if 8 is large enough then the instability cannot develop 
in time and the phase remains essentially real. In the 
opposite limit the fluctuations become large and we have 



Here we will try to examine the role of quantum fluctu- 
ations. Before doing actual calculations, let us give some 
qualitative discussion. As we mentioned before we are 
interested in the regime, where N is large and interac- 
tions are relatively weak A -C N 2 so that the system is 
far from the Mott insulating transition and the quantum 
fluctuations are intrinsically small. This means that nor- 
mally it is possible to use GP approach or at most the 
Bogoliubov extension. However, this is not the case for 
our problem. Indeed, near the classical instability the 
starting point of unstable equilibrium for the Bogoliubov 
expansion of the uniform condensate becomes bad. The 
other way to describe this is to note, that the Bogolui- 
bov equations are nothing but the quantized version of 
linearized equations (|11II12|> . which can predict the onset 
of the instability but fail to describe the nonlinear regime. 
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On the other hand, we can anticipate that the quantum 
fluctuations will remain weak until we cross the instabil- 
ity point. After that they will force the system to evolve 
into the superposition of unstable classical trajectories 
and become unimportant again, when those trajectories 
will be relatively far from each other. 

These ideas known as a truncated Wigner approxi- 
mation^ have been recently applied for the description 
of BECa 18 ! 19 ! 20 ! 21 ! 22 . The usual method of deriving this 
scheme is based on the cubic Fokker-Planck equations of 
motion for the density matrix written in the Wigner rep- 
resentation. The main difficulty with the these Fokker- 
Planck equations is that the third order derivative terms, 
which are responsible for the corrections to GP trajecto- 
ries, makes them practically intractable^. Another en- 
tirely different Keldysh technique can be also used to 
attack non-equilibrium problems 28 . Classical equations 
of motion appear there as the saddle point of the ac- 
tion. However any attempt to include quantum fluctua- 
tions immediately results in a path integral with a self- 
interacting classical field 2 ^, which is very complicated. 
In the Appendix we will show how the TWA nat- 
urally arises from the path integral formulation of the 
dynamics and emphasize the key difference between the 
present derivation and that of the conventional Keldysh 
technique. Moreover, in Ref.[27j we show that it is 
straightforward to go beyond TWA perturbatively in- 
cluding quantum effect on the classical trajectories them- 
selves. We will also argue there that the TWA gives exact 
short-time asymptotical behavior of the evolution of any 
system, the time when it breaks down, however, depends 
on the details of the particular process (see Sec. IIII B|) . 

The whole idea of the TWA is that the expectation 
value of any given operator SI at time t is equal to the cor- 
responding classical observable Sl c i(t) evaluated accord- 
ing to standard GP equations and averaged over an en- 
semble of initial conditions distributed according to the 
Wigner transform of the initial density matrix (see Ap- 
pendix^ for the details of the derivation): 

Sl(t) = J d^d^p^^Sl^t),^ (*)»*). ( 25 ) 
where p is defined as: 

PW>o,V>o) = J dr]odrj {ijjo - ylpolV'o + y) 

x e -Wol 2 -iM 2 e h(no->Po-vo4>o) _ (26) 

In the equation above \ipo±r]o/2\ denote coherent states. 
We use the following measure 



d^ttpQ a d^SlpO a 



(27) 



with the product taken over continuous or discrete spa- 
tial indices, which we suppressed in and (|2l))l to 
shorten the notations. The interpretation of p(V'O)^'o) 
as a probability is not very precise, because the Wigner 



transform does not have to be positive. To get the func- 
tion fi c ((V>> ip*, t) we need to rewrite the quantum opera- 
tor St in the fully symmetrized form and substitute field 
operators a and a' by their classical counterparts ip and 
ip*. In particular, the relation between Sl c i and a more 
familiar version of the classical counterpart of the normal 
ordered operator Si usually appearing in the functional 
integrals is: 



Sl cl = (Sl(ip* + r]*/2^-ri/2)), 



(28) 



where the average is taken over 7/ with the weight 
exp(-|r?| 2 /2). 

Let us now give general comments on validity of l|25fl . 
If the Hamiltonian is non-interacting, then this expres- 
sion is exact. So TWA includes the Bogoliubov approx- 
imation and goes beyond. To recover the latter we just 
need to linearize the classical GP equations of motion 
while evaluating ip(t). This statement is not surprising 
since the noninteracting evolution is always identical to 
classical2ii22*^. If there are nonlinear interactions, then 
in general, there will be corrections to the equations of 
motion themselves. We consider them in Ref . pT;| . Let 
us only note here that for the two coupled condensates 
we showed in Ref.fl^ (see also Sec. 1111 B*|) that the time 
where GP breaks down in the worst possible scenario 
with the least classical initial state having completely 
undefined phase is equal to t c w N/X = J/U. We ex- 
pect that the scaling with N is generic^ and therefore 
(|25|l should be valid at least for the times shorter than 
t c . In next section we will see that if there are no sudden 
perturbations, so that a small fraction of quantum lev- 
els is populated then the time scale of validity of TWA 
becomes much longer. 



Two coupled condensates. Comparison with 
exact solution. 



The main purpose of this section is to test the trun- 
cated Wigner approximation on a simple example of two 
coupled condensates, where it is straightforward to ob- 
tain the exact solution. The two-well version of the 
hamiltonian reads: 

Hit) = -Jalrfaf} + ^ a^a^ - 1), (29) 

where a, j3 = L,R denote the right or the left well respec- 
tively, t" 13 , a = x,y,z are the Pauli matrices. As usually 
we imply implicit summation over repeated indices. Be- 
cause of total number conservation (|29|l is equivalent to 
a more familiar version 

H(t) = -Jairfap + ^ (4 T °%) 2 , (30) 

A convenient choice of the observable is 

n = ^(4r«%) 2 = : {alrfapf : +1 (31) 
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where semicolons denote the normal order. The operator 
n is nothing but the variance of the relative number dis- 
tribution. Let us consider several examples of evolution: 
(i) the initial state is symmetric and the interaction in- 
creases with time, (ii) the initial state is antisymmetric 
and the interaction increases with time, and (hi) the ini- 
tial state is the Fock state and the interaction does not 
change in time. The situation (ii) is directly relevant to 
the "cat" state dynamics we consider in this paper, but 
we also look to the other possibilities to check the validity 
of this approach in a more general case. 

The classical function £l c i(i/j, ip* , t) can be either found 
from the normal-ordered form of the operator Q accord- 
ing to (|28|l or by direct symmetrization of the latter. In 
our particular case it reads: 



1 



N N 8N 2 
1 



87V 2 



(32) 



The final step is to find the probability function p(ipo, ip^) 
according to (|26|l . This will depend on the details of 
the state \0), therefore we have to study different initial 
configurations explicitly. 



1. Symmetric or antisymmetric initial state 

Suppose that at t = the interaction was negligi- 
ble. Then the products of symmetric and antisymmetric 
wavefunctions give the ground and the most excited sta- 
tionary states, respectively. They can be also represented 
as a superposition of products of the two coherent states 
with equal or shifted by tt phases: 



|0) = (4ttA0 



l/4 e -AT 



d6 
— e 

2tt 



-2i6N 



Ne 



id \ 



Ne 



iO-\-i7r<T\ 



(33) 

where a = 0, 1 for the symmetric or antisymmetric state, 
respectively. The integral over the global phase 9 en- 
sures the particle number conservation. Before proceed- 
ing with further analysis let us look into a simpler exam- 
ple of just a product of the two coherent states, where 
the global phase symmetry is broken and 9 takes some 
particular value. Then after straightforward calculation 
one can show that: 



pW>o^o) = 4e- 2 £JV>o Q -^c 



(34) 



We see that in this case the probability distribution of ipo 
is just a gaussian centered near the classical value with 
the relative variance of fluctuations of the order of 1/ \/~N. 
This is completely reasonable and we indeed recover GP 
picture having a single initial state in the limit N — > oo. 
Now let us look closer to the wave function 13311 . After a 
simple calculation the final expression for the probability 
p reads: 



-|Vo + | 2 -|^o-| ; 



where ipo± = V'o l ± V'o R m the symmetric state, and we 
should interchange -00 + and V'o- f° r the it state, L2n(x) 
is the Laguerre's polynomial. This expression is a gaus- 
sian in terms of | Vo ! > however it has a nonlocal behavior 

as a function of IV'o + l- Moreover p^OiV'o) * s n °t posi- 
tively defined. This leads to interesting consequences. 
For example, while the average of IV'o + l 2 , computed with 
help of l|35|) gives the expected classical result, the vari- 
ance of IV'o + l 2 is negative. In Fig. ^ we plot the nor- 
malized function IV'o + l pGV'o + I) f° r the situation with 8 
bosons per well. The extra factor of IV'o + l comes from 
the integral measure: 



IVo + bdV'o + IMVo + l = 1, (36) 



For convenience we rescaled the fields V'o — ¥ vN ip so that 
the classical expectation value of IV'o + l is 2. In the limit 




0.0 0.5 



FIG. 1: Distribution of the GP initial conditions versus IV'o + l 
for the symmetric state with eight noninteracting bosons per 
well. 

TV — > oo we again recover the classical result (for the 
rescaled fields) -0ol,_r = 1 or IV'o + l = 2, 1^0-1 = 0, but 
in a peculiar way. The contributions from | V"o + 1 < 2 will 
cancel each other because of fast oscillations of the prob- 
ability p and only the small interval around IV'o + l = 2 
will give the contribution to the final result. We might 
think, that if the observable is a smooth function of the 
initial parameters, then the details of the distribution 
p(V'0;V'o) are n °t important and we can substitute it by 
some gaussian function with appropriate mean and vari- 
ance. However, as we pointed out before the variance 
given by the distribution is negative, so at least this 
is not very straightforward to do. 

Now let us assume that the interaction increases with 
time according to: 



-^2^(2^0- 



(35) 



A(i) 



tanh(<5t) 
1 - St 



(37) 
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N=8 

0.07-, 




■5 0.02- 



0.01 - 

0.00 -I , 1 , 1 , 1 , 1 , 1 , 

1 2 3 4 5 6 

Interaction X 



FIG. 2: Dependence of the number variance on the interaction 
changing with time according to (13 7t for initial symmetric 
state and 8 bosons per well. Dashed and solid lines show semi- 
classical and exact solutions, respectively. The lower graphs 
correspond to a faster turning on the interaction. 

where S -C 1 is the adiabatic parameter. This depen- 
dence is somewhat different from what we used in the pre- 
vious section. But the resulting instability is still there, 
and besides the main purpose of this section is to test our 
approximation scheme rather then to do some particular 
calculations. The resulting graphs for both symmetric 
and antisymmetric initial states are plotted in Figs. [21 
and 03 Note that even for the eight particles per well 
the agreement between the exact and the TWA solutions 
is remarkable. For the thirty two particles there is a 
small discrepancy for the intermediate value S. Appar- 
ently the semiclassical curve does not capture the small 
oscillations very well. But note that both in the limit 
of large and small 8 the oscillations disappear and the 
agreement becomes perfect. 

Notice that the steady state for the initial antisym- 
metric conditions is exactly the maximally entangled 
"Schrodinger cat" state, where all the bosons occupy ei- 
ther left or right well: 

\9 f ) = -j=(\LLL...) + \RRR...)). (38) 

The ultimate reason for this is that as we mentioned 
above the n shifted state is at classical equilibrium for 
any interactions A. This equilibrium though becomes 
unstable for A > A c . So any fluctuation will cause a 
classical trajectory to end up either in the left or in the 
right well and the quantum zero point motion gives us 
these fluctuations. However, quantum mechanically we 
do not break the left-right symmetry, because it is the 
property of the full hamiltonian. The only way to rec- 
oncile these two results is to have the final configuration 
in the coherent superposition of the left and the right 



a) N=8 



1.0- 




o.o -I — , — , — , — , — , — , — , — , — , — , — , — I 

1 2 3 4 5 6 

Interaction X 



b) N=32 

.0-t— ^ 




Interaction X 



FIG. 3: Same as in Fig.|5|but for initial antisymmetric state. 
The graphs (a) and (b) correspond to 8 and 32 bosons per 
well, respectively . The inset on the graph (b) corresponds to 
the slowest evolution. 



states. This statement can be also verified numerically. 



2. Initial number state 

Here we revisit our results derived earliest assuming 
the two condensates are initially uncoupled and their 
wavefunction is just a product of the two number (Fock) 
states. Then at t = the tunnelling is suddenly turned 
on and the number variance starts to experience some 
oscillatory behavior—. Repeating the same analysis as in 
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the previous section we find: 

/•CfO />OC />2lT 

((4(i)r z Q <%(i)) 2 ) = / / / dn L dn R d9 
Jo Jo Jo 

Pnum(n L )p (n fl )(^*(t)7?'%(t)) 2 , (39) 

where til,r — \ipL,R.(t = 0)| 2 , 9 is the initial phase dif- 
ference between ipr, and ipR. So in the Fock state the 
phases in the two wells are indeed uncorrelated as we 
argued in Ref.|12|. However the number of bosons is dis- 
tributed according to p n um (n) given below and not fixed 
at n = N as we might naively think. In (|39fl we ignored 
an additive l/8iV 2 correction (see i[32|l). The probability 
of having initial occupation n in either well is2&: 

Pnum(n) = 2e- 2n L N (4n), (40) 

where as before Ln(x) stands for Laguerre's polynomial 
of the order N . The function p nU m is very similar to 
its counterpart defined in (|35JI in the sense that it has 
also an oscillatory behavior for n < N and exponen- 
tially decays for n > N. In Ref.^l we showed that 
the simple GP picture (i) gives a multiplicative error 
(1 + 2/N) in the number variance even in the noninter- 
acting limit, and (ii) it is valid for the finite amount of 
time shorter than some characteristic scale determined 
by interactions: t < t c w J/U = N/X. For longer times 
the GP result starts to deviate strongly from the exact 
solution due to recurrence occurring in a quantum sys- 
tem. We might guess that the agreement between the 
semiclassical and the quantum results can be improved 
upon including quantum fluctuations at initial time ac- 
cording to l|40l) . This is indeed the case for t < t c , i.e. 
the discrepancy in the prefactor completely disappears. 
However, as we can see from Fig.0J these fluctuations do 
not affect the time t c itself so that the correct result can 
be recovered only if we also include quantum scattering, 
or in other words deviations of the trajectories from the 
classical ones. Finally we would like to note that the ini- 
tial number state is the worst possible from the classical 
point of view, because the phase is completely undefined 
there. So we expect that in general t c gives the lower 
boundary of the applicability of TWA. 

IV. NUMERICAL RESULTS FOR MULTIPLE 
WELLS 

Having established the general framework and checked 
its validity let us move on to the main subject of the 
paper. First, following the analysis given in Sec. [HI we 
will study the temporal behavior in a periodic array of 
wells being initially in the ir state. Then we will consider 
a case of a harmonic trapping potential. 

3. Periodic Array 



N=16 

0.5 -, 




2 4 6 8 10 
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FIG. 4: Time dependence of the number variance for the 
initial number state and 16 bosons per well. The interaction 
strength is A = 1. Note that here and elsewhere in this paper 
time is measured in units of h/ J. 

N bosons per well (we assume M to be even) is: 
|0) = (inNM^e^ f ^- e ~ ieNM ]J l^Ne"*"^, 

J ^ 3 = 1 

(41) 

where | . . . )j stands for the coherent state in the j-th well. 
This is an eigenstate of the noninteracting hamiltonian 
and apart from the global phase 9, which conserves the 
total number of bosons, it is just a product of coherent 
states with alternating phases. Ignoring the integral over 
9, results in a gaussian probability distribution of the 
initial state (compare with Q34[0: 

M _ 
p(4>o,%) = 2 M JJe- 2 ^-^ 4 *'! 2 . (42) 

3=1 

while the correct result for 141|) reads: 

p(i> ,r ) = 2 A V 2M ^l^l 2 W(4M|^ 00 | 2 ), (43) 
where i/jq k stands for the discrete Fourier transform of 

i^of 

^0k = ^J2^je- ikj (44) 

3 

Clearly as M — > oo the difference between (|42H and l)43|) 
should vanish. Next let us define the operator 0, which 
would be a good measure of the instability: 



The straightforward generalization of ll-i-'il) for the 7r 
state in a periodic chain of M coupled condensates with 



n= ^M(M-l) B°^-^ 2 - (45) 
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This is just a normalized sum of number variances over 
the different wells. We have chosen the prefactor so that 
(fi) < 1, with 1 corresponding to the state with all the 
bosons located in any single well. It is easy to verify that 
the classical counterpart of O is: 



1 



N 2 M{M - 1) 
1 

47V 2 (A/ - 1)' 



E m 



N - 



(46) 



It is reasonable to expect that as the number of wells in- 
creases the global phase becomes less and less important 
and therefore (|42[1 becomes more and more accurate. 

Next we consider several specific examples. First let 
us take the interaction to be monotonically increasing in 
time according to Q37|). In Fig. JSJ we plot the result- 
ing evolution of the state for the case of ten wells with 
eight bosons per well. The solid and the dashed lines 
correspond to the probability distributions given by l|4rt|l 
and l|42[l . respectively . Clearly there is no significant 
difference between them. Note that the upper curves 
corresponding to smaller <5, i.e. to the adiabatic limit, 
saturate at Q « 1, which shows that all the bosons ap- 
pear in a single well, i.e. the resulting steady state \F) 
is: 



\F) 



1 



<M 



1111 



1222.. 



\MMM ...)). (47) 



The whole procedure of driving the system into the maxi- 
mally entangled state described here is co ncep tually very 
similar to that recently suggested in Ref.^llll, where the 
tunneling was assumed to decrease by the spatial drag 
of the double-well condensate through a beam splitter. 
The important difference however, is that here we are 
not limited by the double-well system and can consider 
larger arrays, so that our entangled "cat" occupies more 
than two macroscopic states. 

There is an important issue, which was completely ob- 
scured in the preceding analysis. Indeed, studying the 
number variance alone, it is impossible to distinguish the 
"cat" state from the collapsed condensate. While the 
collapse is often very well reproduced using GP equa- 
tions, it is much harder to describe the recovery within 
this framework. To examine this issue let us consider the 
interaction, which is periodic in time: 



A(i) = A sin 2 (7r^), 



(48) 



where the parameter S as in (|57|) det ermines the adia- 
baticity of the process. If S is small then we expect com- 
plete restoration of the initial state after one period of 
oscillation T = 1/5. With decreasing period we gradu- 
ally loose adiabatic limit and the evolution of the system 
is no longer expected to be periodic. Fig. summarizes 
this discussion and the graphs are in perfect agreement 
with our expectations. The phase correlation in this fig- 
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FIG. 5: Number variance as a function of interaction for the 
case of ten wells with eight bosons per well. The solid and the 
dashed lines correspond to the distributions given by (I43H and 
(1 121 . respectively. The upper curves correspond to a slower 
increase in interaction. 



ure is defined in a usual way as 
1 



2NM 



^2(a\a 3+1 ) +c.c 



and the "— " sign is inserted for the sake of convenience 
because we start from the 7r-state. 
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FIG. 6: Number variance and the phase correlation as a func- 
tion of time divided by the period of interaction (T = 1/5) for 
ten wells and eight bosons per well. The interaction changes 
with time according to 1481 . For smaller 8 the phase restora- 
tion is almost complete, which proves the coherence of the 
dynamics. 

So far we considered examples, where the initial state 
was a noninteracting 7r-phaseshifted condensate and then 
the interaction was slowly ramped up (or equivalently 
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tunnelling was ramped down). This rather hypothetical 
situation is certainly suitable for the theoretical analysis, 
but hardly applicable to any real experimental realiza- 
tion, where the interactions between the bosons is al- 
ways present. So the initial wavefunction is not a simple 
product of coherent states but a more complicated ob- 
ject. One way to proceed with this issue is to write down 
an approximate wavefunction, which can be obtained us- 
ing either Bogoliubov hamiltonian (valid for A <C iV 2 ) or 
variational methods, and then find its Wigner transform 
according to J2HJ). This would be a tractable but lengthy 
calculation. Instead we can do a simple trick, demon- 
strating the advantage of the present approach. Namely, 
we can start from the simple noninteracting ground state, 
then adiabatically increase the interaction strength to the 
desired value, and finally apply the 7r-phaseshift and fol- 
low the subsequent evolution. In this scheme there is no 
need to do any additional analytic calculations, both the 
"fake" evolution to the interacting ground state and the 
subsequent real dynamics are described within the same 
scheme. Moreover the computational time does not in- 
crease much because of the extra "fake" evolution to the 
true ground state. Indeed, the classical motion is stable 
before the 7r-phaseshift is applied so that GP equations 
can be efficiently integrated. To be more specific, as- 
sume that we start from the interacting condensate with 
A = const(t), then at t = suddenly apply a n phase- 
shift and follow the subsequent appearance of the "cat" . 
Clearly if A > A c (see (JSJJ) the system will become un- 
stable and the quantum fluctuations will force it into the 
entangled state. Without any calculations, it is obvious 
that the final state will not be maximally entangled as 
described in Fig. because of the energy conservation. 
Indeed, the state with all the bosons occupying one well 
costs a huge interaction energy which can not be com- 
pensated unless there is an external pumping resulting in 
the time dependence of the hamiltonian. From the same 
considerations it is obvious that the maximum possible 
entanglement within this scheme can be achieved at some 
intermediate values of A. Thus if A < A c , there is no in- 
stability to begin with and so no entanglement in the final 
state. On the other hand for the large A any significant 
number fluctuations will cost a strong increase in the in- 
teraction energy, which can not be compensated by a lim- 
ited decrease of the hopping energy. We plot the result- 
ing number variance as a function of time for the periodic 
array of ten wells in Fig. [7] The critical value of the inter- 
action for this case is A c = 2sin 2 (7r/10) « 0.19. Clearly 
for A getting closer to A c the oscillations become slower 
and the steady state is achieved later. We would like 
to point out, that this type of dynamics certainly corre- 
sponds to a sudden perturbation discussed in Sec. lIII B"2l 
Therefore we expect a finite ergodic time for a full quan- 
tum solution so that there is no true steady state for the 
finite number of bosons. Therefore to achieve a stable 
stationary entangled state it is always preferable to drive 
the system adiabatically towards the instability, i.e. ap- 
ply a 7r-phaseshift at smallest possible A and then slowly 



ramp up the interaction. 




FIG. 7: Number variance as a function of time for ten wells 
and eight bosons per well. The interaction remains constant 
in time in this example but at t = a sudden 7r-phaseshift is 
applied to the system. 



4- Harmonic 



The effect of harmonic trap can be mimicked by adding 
a quadratic potential term to the hamiltonian of the sys- 
tem |JTJ: 



v, = 



(49) 



Here ujt is the trapping frequency. 

As before let us assume that initially the interactions 
are suppressed and the system is in the ground state: 



N tot 



\o)= 5>*°S 



I Vac), 



(50) 



where |Vac) denotes a vacuum state with no particles, 
Ntot is the total number of bosons in the system and 



V2 

(iTLJt)' 



(51) 



is the ground state wave function of a single boson in 
the harmonic trap in the coordinate representation. The 
state \0) can be also written as: 

dO 



|0)= / ^ el9{ni+n2+ '''~ Ntot) \^^ a ^^^ta 2 )c..., 

(52) 

where as usual \a) c stands for the coherent state. If the 
number of populated wells is not small, then we can ig- 
nore the global phase and use an approximate expression 
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\0)*H\VNi 



totOtj 



(53) 



As in the previous discussion we assume that at initial 
time the phases in the adjacent wells were uniformly 
shifted by some phase <6. However, we will not consider 
only the case with <f> = tt, because the tt phase shift in 
the nonuniform potential does not give a stationary solu- 
tion, although it still preserves inversion symmetry. Also 
the number variance is no longer a convenient measure 
of the instability, because the distribution is not uniform 
even in the initial state. Instead let us introduce two 
other quantities: coordinate of the center of mass and 
the condensate width: 



-X 2 . (54) 



The semiclassical operators can be trivially obtained 
from l|54(l . For example, 



^ = iEj(^,-l/2)>. 



(55) 



Let us assume that both the interaction and the trap 
frequency increase in time: 



A(i) 



tanh 5t 
1-dt ' 



"?(0) 
l-5t' 



(56) 



We define A according to (0) with N = N tot \ao\ 2 be- 
ing the number of bosons in a central well. Note that 
with u>t there are two degrees of freedom: A/ J and 
Lo t / J. If we want to reflect the experimentally relevant 
case of vanishing tunnelling rather then increasing inter- 
action we have to keep the ratio \/ui 2 constant. Extra 
multiplier "tanh 8t n in the interaction term is taken only 
for convenience purposes and it does not qualitatively 
change any of the results. 

The two graphs in Fig. [H] show the condensate width 
and the center of mass position versus time for differ- 
ent initial phase imprints. Note that if the phaseshift is 
equal to tt, there is no center of mass oscillation because 
of the inversion symmetry. However the width oscilla- 
tions are very large and pronounced in this case. They 
also die very fast away from the tt phaseshift, however 
the steady state condensate width depends on <f> rather 
smoothly. It is also interesting to note that the center 
of mass oscillations decay faster for larger angles. This 
is in the direct agreement with the previous predictions 
of self-trapping at large phase-gradients^. We would like 
to point out that the trapping is real in our case because 
of the rapidly increasing interaction (or equivalently van- 
ishing tunnelling). 

There is, however, a major difference between the re- 
sults for a periodic array and a harmonic trap. In the 
classical limit the 7r-state is an eigenstate in a periodic 





FIG. 8: Center of mass position and the condensate width as 
a function of time for the interaction and trapping frequency 
increasing in time according to 156|l with 8 = 0.01, uit(fl) = 
0.03, and the total number of bosons N to t = 100. 



array for any strength of the interaction, therefore the 
motion is completely driven by the quantum-mechanical 
fluctuations, while in a harmonic trap the tt state is classi- 
cally unstable (that is why we plot our observables versus 
interaction for the periodic array and versus time for the 
parabolic trap). The classical instability makes quantum 
effects unimportant if the number of bosons is large. To 
further elucidate this point let us compare the results for 
a single Gross-Pitaevskii trajectory, i.e. simple classical 
limit, with the TWA result. The two dependences are 
shown in Fig. [5] and they are clearly very alike. This 
proves that the resulting instability has a simple classi- 
cal nature and does not give a "cat" state, rather the 
condensate simply has large amplitude breathing modes. 
Such an outcome should not be very surprising. In order 
to get to the macroscopical entanglement one has to pre- 
pare a classically equilibrium but unstable state, which 
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Time 



FIG. 9: Condensate width as a function of time for the same 
conditions as in Fig. @. The dashed line is the simple Gross- 
Pitaevskii result and the solid line is the improved calculation 
according to l|25ll . 
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FIG. 10: Phase coherence as a function of time for the pe- 
riodic interaction I48H with S — 0.004. The total number of 
bosons Ntot = 100, the trap frequency u>t = 0.03, and the 
initial phaseshift cj> = n. Clearly the dynamics is completely 
irreversible as opposed to Fig.HJ 



can be driven away via quantum fluctuations. This can 
be achieved by a simple phase shift only in a uniform 
periodic lattice. 

We can also check the reversibility of the evolution for 
the periodic in time interaction (|48|l . And clearly, (see 
Fig. Ill) [1 the dynamics is completely irreversible. This re- 
sult is also not surprising given that the initial state is not 
stationary even without interactions, so the 7r phaseshift 
looks as a sudden perturbation applied at t = and even 
if the interactions change with time slowly the system is 



not in the adiabatic regime. 

V. DISCUSSION 

In the previous sections we adopted a reduced Wigner 
approximation allowing to treat dynamical evolution of 
classical instabilities due to quantum fluctuations. In 
this analysis we completely ignored the classical noise, 
which may come from various sources e.g. fluctuations 
of the laser intensity or wavelength, fluctuations of the 
magnetic field creating the harmonic trap, collisions with 
external atoms, triple interactions between the bosons 
leading to the loss of the atoms from the condensate, etc. 
All these sources are generally quite weak, otherwise no 
beautiful experiments would ever be possible. On the 
other hand the "cat" states are extremely fragile to any 
kinds of classical noise. For example, we know that in a 
double well in the equilibrium the noise required to de- 
stroy the coherence exponentially scales with the mass 
(or equivalently number of particles). So that heavy ob- 
jects are always localized in one of the classical minima. 
However, this is not the case in our situation. Indeed, the 
system is far from the equilibrium and the "cat" state is 
not a ground state. Let us crudely estimate the limits on 
the noise. Assume that the classical fluctuations result 
in a random external potential with a usual correlator 
(Y(t)Y(t')) — Yo5(t — t') and consider the case of two 
wells for simplicity. Then the relative phase flow with 
time due to fluctuations will be 

59 ~ NYoVt- 

On the other hand the minimum time required to get to 
the cat state is (see Sec. nUl : 

t ~ - ~ In A. 
o 

Requiring that the total phase accumulated during the 
evolution is less then one we get the upper bound for the 
noise which does not destroy the coherence : 



AVhiiV' 

This result is encouraging, because instead of exponen- 
tional we get a much weaker scaling of the noise with 
the number of particles. Clearly, in the multi-well case 
there will be an extra scaling with the number of wells 
and of course no cat state is possible in the thermody- 
namic limit. However, the scaling will be again weak and 
tractable. Another possibility to observe macroscopic 
"cat" states experimentally is to use a modulated hop- 
ping between the sites in a large array effectively split- 
ting the condensate into pairs of sites, so the effects of 
the classical noise become weaker. 

The other two constraints used in our analysis that 
the number of bosons is large and the condensate is one- 
dimensional are also not essential. We used the large 



12 



number of bosons rather to satisfy the formalism then 
to explain the effect. The instability is there even in 
the fully quantum-mechanical treatment of the problem 
and the resulting maximally entangled state is clearly in- 
dependent of N. The only thing is that the evolution 
of macroscopically entangled states with a large num- 
ber of atoms is more interesting from the experimental 
prospect. Concerning higher dimensions we would like 
to point out that, although we did not performed actual 
calculations, it is extremely unlikely that any of the re- 
sults will qualitatively change. The instability will clearly 
survive in any dimensions and the 7r-state can be imple- 
mented, at least theoretically, in square lattices of arbi- 
trary dimensionality. 

In conclusion, we showed that it is possible to cre- 
ate macroscopically entangled "Schrodinger cat" states 
in the bosonic systems in optical lattices with finite num- 
ber of sites. We justified and used truncated Wigner ap- 
proximation generalizing a simple GP approach by the 
exact treatment of quantum fluctuations at initial time 
of the evolution. The resulting expressions were tested 
on a solvable model of the two coupled condensates and 
we found a very good agreement with the exact results. 
At the end we presented some numerical simulations for 
the multiple-well condensates and argued that it is pos- 
sible to create a "cat" state with more than two possible 
outcomes. 
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APPENDIX A: MICROSCOPIC DERIVATION OF 
THE TRUNCATED WIGNER APPROXIMATION 

Let us assume that at the initial time of evolution t = 
the system is described by some density matrix po : 

Po = £ P et)lxXxl, (Al) 

X 

where \\) represents some basis and P(x) is the proba- 
bility to be in a particular state (if we use coherent ba- 
sis then P(x) coincides with Glauber P-function widely 
used in quantum optics^). If the initial state is pure, 
the sum contains only single term. According to stan- 
dard quantum-mechanical formulas the expectation value 
of an arbitrary normal-ordered operator f2 at arbitrary 
time t is given by: 

fl(t) = Tr [ Po T KT e l ti H ^ dT ^-^o j , (A2) 

where Tk t is the time-ordering operator along the 
Keldysh contour going from to t and then returning 
back to 0, Ti is the Hamiltonian of the system, e.g. in 
our case. In the same way one can define correlation func- 
tions of products of operators. The conventional trick to 
deal with expressions like (|A2I) is to rewrite them in the 
path-integral form using the coherent-state representa- 
tion^. The only difference with the more usual equilib- 
rium case is that there are two exponents containing TL. 
So it is convenient to introduce two fields a / and ab prop- 
agating forward and backward in time 2 ^. Then instead 
of (|A2|) we can write 



Q(t) = J Da f Da b (a b0 \p \af Q ) e -a}o«fo+*} *fi+iH( a } , afl )AT _ e -a* JQ a fl 
fl(a*fQ, QbQ, t)e a *f Q abQ 'e~ a * b Q ab « + °& Q a bQ-i- in ( a t a a t >q-i)At - _ _ e -< »'» i 



(A3) 



where At = t/Q and Q — > oo, Q(a*^ ab,t) is the normal 
ordered operator Q with the fields aJ(t) and a(t) sub- 
stituted by the complex numbers aj(t) and a,b(t). The 
expression above is intentionally written in the discrete 
form, since we want to take special care of the bound- 
ary effects. In the classical limit the evolution equa- 
tions are deterministic, therefore a/ c ;(r) = ab c i{r). On 
the other hand, because of quantum fluctuations the 



two trajectories may be different. Therefore instead of 
the forward and backward fields it is natural to intro- 
duce their classical (ip) arL d quantum (j]) combinations 28 : 
a.f = ip + ri/2, ab = ip~ 7 ?/2- So that in the classical limit 
i/j should correspond to the solution of the GP equations 
and r\ should be simply equal to zero. Now we can take 
a continuum limit in i|A3|l taking At — > so that: 
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e -|V'o| 2 -i| J )o| 2 -ihWI%^(%>o-r,o^) e / t dr t; *(T)£[V.^*,r]-, ) (r)£*W,'A%T] dr ^ (^-* (r )r,(t)+V(r)r,* (r)) |r,(r ) 



where £["0, ip*, t] stands for the classical (GP) differential 
operator acting on the field ip(t): 

C 3 [0, r, T \=^ + (V-i+i + V^-i) - A(r) (A5) 

Note that £ as well as fields ip and 7/ contain spatial in- 
dices which we suppressed in (|A4ll to simplify notations. 
A closer look to equation l|A4(l shows that there are lin- 
ear, quadratic and cubic terms in the quantum field rj the 
latter appearing only due to interactions. It is intuitively 
clear that in the classical limit tj(t) should be small in 
some sense. Thus if we ignore completely all nonlinear 
terms in rj, then the functional integral over the quan- 
tum field becomes a trivial product of <5-functions en- 
forcing GP equations on the classical field tp. The next 
approximation will be to leave quadratic corrections but 
ignore cubic. From (| A4|> it is clear that the quadratic 
corrections affect only the initial and the final times of 
the evolution, that is why it was important for us to start 
from a discrete version of the path integral and be careful 
about boundaries. The integral over rj(t) transforms the 
operator f2 into fl c i according to (j28(l . It is a simple ex- 
ercise to check that 51 c ; is obtained by first symmetrizing 
f2 and then substituting the operators a and by the 
c-numbers ip and ip*. For example if 

ft = a} a = i(a f a + aa f ) - i (A6) 

then 

n cl = (r + rf/2, - 77/2) =n>~\ (A7) 

where as it follows from (|A4|I , the average over r\ is taken 
with the weight exp(— |t7| 2 /2). 

The second quadratic contribution originates from the 
field r/o corresponding to the initial time of evolution. Be- 
cause of the coupling to (^o) in 1|A4(1 . these fluctuations 
introduce a probability distribution for the classical ini- 
tial conditions given by J2HJ). If we ignore the corrections 
to the classical equations of motion coming from the third 
power of the quantum field 77, then the time dependence 
of the observable will be given by i|25|l . 

Let us now give general comments of validity of i|25|) 
and i|26|) . If the Hamiltonian is non-interacting, then 
these expressions are exact. If there are nonlinear inter- 
actions, then in general, there will be corrections to the 
action involving terms proportional to all odd powers of 
77: rf in our case (see or higher in general. Those 

terms will affect the time evolution, which will not be de- 
scribed by the GP equations any longer. In Ref.[2]J we 



show that the corrections to the GP dynamics arise in the 
form of the quantum scattering events, which are equiva- 
lent to the nonlinear response to the infinitesimal pertur- 
bation of the field ip along its classical path. We only note 
here that these corrections are always of the form f(t) /N 
with f(t) being some time dependent function satisfying 
f(t) — > at t — > 0, and 1/N is our semiclassical parame- 
ter. So the TWA given by and always gives the 
exact short time asymptotical behavior of the evolution. 
As we discuss in Sec. IIII Bl the time when TWA breaks 
down depends on the details of a particular process and 
becomes longer under slow perturbations, where only a 
limited number of quantum levels are excited. 

Another point we would like to make is that though the 
expansion in powers of 1/N is clearly around the classi- 
cal solution, there is no h present anywhere. This should 
not be surprising since here and quite often in the atomic 
physics the Planck's constant either completely absorbed 
into energies, which are measured in Hz, or into time. In 
conventional units ftr 1 appears as the prefactor in the 
action justifying the saddle-point or classical approxima- 
tion. In the same way the number of bosons per site 
N appears as a prefactor in the exponent of (|A4(1 after 
rescaling ip — > \fNip and r\ — > \[Nr\. So in general any 
expansion in powers of 77 is in fact the expansion in pow- 
ers of h. 

Let us finally spend a few words discussing the dif- 
ference between the present derivation and that of the 
conventional Keldysh technique. The key point in our 
discussion is that we ascribed the time dependence to 
the operator f2 itself, while leaving the density matrix 
time independent. This allowed us to completely sep- 
arate initial quantum fluctuations, which entered in the 
form of Wigner distribution of initial conditions to classi- 
cal trajectories (|26|l . from the quantum dynamical effects 
(which we consider in Ref.|27jL On the other hand in the 
Keldysh technique the density matrix acquires time de- 
pendence and the initial density matrix is absorbed into 
the quantum propagator^. While it is still possible to 
derive GP equations in the saddle-point approximation, 
integrating out quantum fields in the lowest order gives 
a complicated self-interacting classical action^, which is 
hardly possible to deal with except perturbatively or us- 
ing stochastic methods. This should not be surprising 
since any diagrammatic technique uses a noninteracting 
Gaussian limit as a starting point. Therefore to get just a 
classical GP dynamics in the Keldysh technique, it is nec- 
essary to sum all diagrams with classical vertices (three 
classical fields and one quantum)2&, which looks virtually 
impossible. 
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